import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.optimize import curve_fit

# 1. DATOS DE ENTRADA
v_inf = np.arange(0, 19)
v_sup = np.arange(1, 20)
v_media_int = (v_inf + v_sup) / 2  # Marca de clase (v_i)
horas = np.array([240.8, 682.9, 1019.1, 1208.1, 1243.1, 1155.5, 975.3, 763.4, 
                  555.1, 376.5, 239.0, 141.8, 79.7, 42.0, 21.0, 9.6, 4.4, 1.8, 0.9])

rho = 1.2  # Densidad del aire (kg/m3)
total_horas = horas.sum()
frec_relativa = horas / total_horas
frec_acumulada = np.cumsum(frec_relativa)

# 2. CÁLCULOS DIRECTOS (Inciso a, b, c)
v_media = np.sum(v_media_int * frec_relativa)
v_cubo_media = np.sum((v_media_int**3) * frec_relativa)
v_eficaz = v_cubo_media**(1/3)
densidad_potencia = 0.5 * rho * v_cubo_media

# 3. AJUSTE DE WEIBULL (Inciso d)
# Función de densidad de probabilidad de Weibull
def weibull_pdf(v, k, c):
    return (k/c) * (v/c)**(k-1) * np.exp(-(v/c)**k)

# Ajustamos los parámetros k y c a los datos observados
popt, _ = curve_fit(weibull_pdf, v_media_int, frec_relativa)
k, c = popt

# Parámetros derivados de Weibull
v_media_w = c * gamma(1 + 1/k)
v_cubo_media_w = (c**3) * gamma(1 + 3/k)
v_mediana_w = c * (np.log(2))**(1/k)
v_moda_w = c * ((k-1)/k)**(1/k) if k > 1 else 0
v_eficaz_w = v_cubo_media_w**(1/3)
densidad_potencia_w = 0.5 * rho * v_cubo_media_w
v_cubo_media_w = (c**3) * gamma(1 + 3/k)
v_eficaz_w = v_cubo_media_w**(1/3)

# --- RESULTADOS POR CONSOLA ---
print(f"--- RESULTADOS DIRECTOS ---")
print(f"Velocidad media: {v_media:.2f} m/s")
print(f"Media de v^3: {v_cubo_media:.2f} m^3/s^3")
print(f"Velocidad eficaz: {v_eficaz:.2f} m/s")
print(f"Densidad de potencia media: {densidad_potencia:.2f} W/m^2")
print(f"\n--- PARÁMETROS WEIBULL ---")
print(f"Factor de forma (k): {k:.2f}, Factor de escala (lambda): {c:.2f}")
print(f"V_media (W): {v_media_w:.2f} m/s\nV_mediana (W): {v_mediana_w:.2f} m/s\nV_moda (W): {v_moda_w:.2f} m/s")
print(f"Media de Velocidades al Cubo: {v_cubo_media_w:.2f} m³/s³")
print(f"Velocidad Eficaz (V_rms): {v_eficaz_w:.2f} m/s")
print(f"Densidad potencia (W): {densidad_potencia_w:.2f} W/m^2")

# 4. GRÁFICOS
fig, axs = plt.subplots(2, 2, figsize=(14, 10))

# a) Histograma Frecuencias Relativas
axs[0, 0].bar(v_media_int, frec_relativa, width=1, edgecolor='black', alpha=0.6, label='Observado')
axs[0, 0].set_title('Histograma de Frecuencias Relativas')
axs[0, 0].set_xlabel('v (m/s)')

# b) Frecuencias Acumuladas
axs[0, 1].step(v_sup, frec_acumulada, where='post', color='red')
axs[0, 1].set_title('Histograma de Frecuencias Acumuladas')
axs[0, 1].grid(True, alpha=0.3)

# c) Curva de Duración (v vs horas que se supera esa v)
horas_excedidas = (1 - np.insert(frec_acumulada, 0, 0))[:-1] * total_horas
axs[1, 0].plot(horas_excedidas, v_inf, color='green', lw=2)
axs[1, 0].set_title('Curva de Duración de Velocidad')
axs[1, 0].set_xlabel('Horas al año')
axs[1, 0].set_ylabel('v (m/s)')

# d) Función de Weibull
v_eje = np.linspace(0, 20, 100)
axs[1, 1].bar(v_media_int, frec_relativa, alpha=0.3, label='Datos')
axs[1, 1].plot(v_eje, weibull_pdf(v_eje, k, c), 'r-', lw=2, label=f'Weibull (k={k:.2f}, lambda={c:.2f})')
axs[1, 1].set_title('Ajuste de la Distribución de Weibull')
axs[1, 1].legend()

plt.tight_layout()
plt.show()